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Abstract 

We study stability, dispersion and dissipation properties of four numerical schemes (Iterative 
Crank-Nicolson, 3'rd and 4'th order Runge-Kutta and Courant-Fredrichs-Levy Non-linear). By use 
of a Von Neumann analysis we study the schemes applied to a scalar linear wave equation as well as 
a scalar non-linear wave equation with a type of non-linearity present in GR-equations. Numerical 
testing is done to verify analytic results. We find that the method of lines (MOL) schemes are the 
most dispersive and dissipative schemes. The Courant-Fredrichs-Levy Non-linear (CFLN) scheme 
is most accurate and least dispersive and dissipative, but the absence of dissipation at Nyquist 
frequency, if fact, puts it at a disadvantage in numerical simulation. Overall, the 4'th order Runge- 
Kutta scheme, which has the least amount of dissipation among the MOL schemes, seems to be the 
most suitable compromise between the overall accuracy and damping at short wavelengths. 
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1 Introduction 

The area of numerical relativity is very complex with many factors influencing the calculations. 
Due to this complexity, choosing a best numerical scheme for solving problems in general relativ- 
ity (GR) may be of particular importance in order to have a long term stability and to minimize 
truncation errors. At the same time, problems in numerical relativity usually demand a large 
amount of computational resources. Hence it is also desirable to choose a scheme which best 
exploits the available computer power at any given time, i.e. provides accurate results with a 
minimal number of operations. 

First and foremost the scheme needs to be numerically stable, but in addition to this ba- 
sic (but not always trivially established) property, there are a number of other features, such 
as accuracy, dissipation and dispersion properties, which may have important impacts on the 
numerical solutions as well. 

The purpose of this paper is to investigate numerical properties (in particular dispersion, sta- 
bility and dissipation) of some numerical schemes which can be used to solve GR type equations, 
allowing users of these schemes to estimate potential pitfalls and limitations of the schemes. 
The schemes that are being analyzed are the iterative Crank- Nicolson scheme (ICN), the third 
and fourth order Runge-Kutta schemes (RK3 and RK4) and a nonlinear version of the classical 
Courant-Friedrichs-Levy scheme (CFLN). A unique aspect of numerical GR is that the inte- 
gration may be affected not only by stability of the scheme, but also by constraint and gauge 
instabilities which are intrinsic to the equations themselves. In order to separate these effects 
from the stability of the schemes and to investigate strengths and weaknesses of the different 
schemes, we want to apply the schemes to simpler scalar wave equations. Numerical schemes 
are usually analyzed by their application to the standard linear scalar wave equation (e.g. PP, 
j2]) however, they are usually applied to solve nonlinear problems. By use of a classical Von 
Neumann analysis, we analyze properties of the schemes applied to both the standard scalar 
linear wave equation and to a nonlinear scalar wave equation, which is designed to mimic some 
properties of the nonlinear terms in the Einstein equations. By this we hope to expand our 
knowledge of the behavior of the schemes in a regime which is closer to the computational real- 
ity of numerical GR then a sole analysis of the schemes in the fully linear regime. The non-linear 
wave equation, a perturbation analysis of it and the numerical schemes are presented in section 
121 The results from the Von Neumann analysis of the schemes applied to the scalar linear and 
nonlinear wave equations are presented in section El and El respectively. To support the analytic 
results, numerical results are presented in sectional 

2 Formulation of the problem 

In this section we present the numerical schemes which we will analyze in subsequent sections, 
the non-linear scalar wave equation to which we will apply the schemes and the method by which 
we will be analysis the schemes. 
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2.1 A Non-linear Wave Equation 

Consider the following scalar, quasi-linear, hyperbolic partial differential equation of two inde- 
pendent variables, t and x (0): 



d'^g d'^g 1 [ dg^ ^ 



dt"^ dx'^ g \dt 



(2.1) 



or cast into first order form: 



dg_ 

dt 



K 



d^_d^_ {Kf_ ^^-^^ 
dt dx'^ g 

This equation is essentially the standard scalar wave equation with an added non-linear term. 
The equation is interesting because the non-linear term mimics part of the non-linearities present 
in GR equations. This can be seen by recalling that the structure of a Ricci tensor Rab is 
~ ^ 9r + ^ rr, where Christoffel symbols T ~ g~^dg, and g is the metric. Thus Rab can be 
represented as a sum of the terms R ~ J29^^^^9 + J2 9~^i^9y- Equation (j2.ip thus mimics a 
type of non-linearity present in GR equations Rab = 0. Especially equations (|2.2p resembles the 
evolutionary part of GR equation in a standard ADM 3-1-1 form (with zero shift and constant 
lapse). We expect that in order to successfully apply the schemes to GR equations they must 
be applicable equation ()2.1|) . 

Equation 1)2.11) posses a large number of non-trivial solutions 0. For test purposes we use 
two particular analytical solutions to this equation: A spatial constant solution: 

git) = y/Ci + C2-t, {Ci + C2-t> 0)} (2.3) 

and an exponential solution: 



g{x, t) = exp ±VC ■ X ± \ — ■ t ] , {C > 0, ± signs independent)^. (2.4) 




2.2 Perturbation analysis 

We wish to analyze the analytical behavior of small amplitude perturbations of solutions ()2.3p 
and (El- 

Let go{x,t) be a base solution of equation (|2.1|) and let go + g he a perturbed solution, with 
1^1 ^ \9o\- We linearize ()2.H) around go and obtain a linear equation for the perturbations 

9^9 d^9 f ^ dgo\ dl ^ / l_dgoY ^ (2 5) 



dt^ dx"^ \go dt J dt \go dt 

^di't) ~ ^VCi + C2 ■ t is naturally also a solution but since we are trying to mimic some features of the Einstein 
equations, g (which mimics the 3-metric of the general relativity) cannot be negative. This is also why we require 
Ci + C2 ■ t > as we would otherwise obtain complex solutions which would be unphysical. 

■^C cannot be negative as this would create complex solutions (see previous footnote). 
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or 

^-^^2A^ + A'd (26) 

^--t (2.7) 
Consider perturbations of the form: 

g oc e^^'-^^\ (2.8) 
where / = Then the dispersion relation is 

= ^2 ^ 2AIuj - A^. (2.9) 
Substituting u = ur + Iuj into equation ()2.9|) and separating real and imaginary parts we obtain: 

2IujRUJi = 2IAujR {2 AO) 

uji-uj = -2Aui- A^ (2.11) 
From (j2.1(Jj) it follows that uj = A, hence the growth rate of a perturbation g oc e^'^'^'^ is 

g oc e"^*. (2.12) 

We see that a base solution is stable if A > and otherwise it is unstable. Hence growing 
solutions are stable and decaying solutions are unstable. 

Due to a local linearizion involved in the above analysis, the conclusion is strictly valid for 
perturbations with high wave numbers k = ^ compared to the base solutions, and only if the 
growth rate is ^ then that of the base solution. By looking at the spatial constant base solution 
fl2.H|l . perturbed by a spatial constant perturbation {k = 0), we can estimate the behavior of 
spatial constant or very long wavelength perturbations. 

Consider a spatial constant base solution ()2.3|) perturbed by a small spatially constant per- 
turbation: g(t) = go{t) +g{t). Then equation ()2.5|) simplifies to an ordinary differential equation 
which can be solved to yield: 

where K is an arbitrary integration constant. From ()2.13|) we see that for growing base solutions 
(i.e positive C2 in ()2.3p ) the second term is growing proportionally to the base solution, while 
the first term is decaying inversely proportionally to the base solution, hence the relative error 
is decreasing towards a constant value proportional to K. That is, the perturbation is growing 
but it is not unstable. If the base solution is decaying (i.e. negative C2 in ()2.3|l ). the first term 
increases rapidly as go{t) —>■ and the solution becomes unstable as the relative error grows 
unbounded. 



where 
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2.3 The Numerical Schemes 

The four schemes that we are investigating are 

• Iterative Crank-Nicolson 

• 3'rd order Runge-Kutta 

• 4'th order Runge-Kutta 

• Courant-Friedrichs-Levy Nonlinear 

The first three schemes are based on the method of lines approach P], while the Courant- 
Friedrichs-Levy Nonlinear scheme is based on the classic central-difference, second-order explicit 
scheme introduced in 1928 by Courant, Friedrichs and Levy [S]. The schemes are defined as 
follows: 



2.3.1 The Iterative Crank-Nicolson scheme 

The iterative Crank-Nicolson scheme (ICN) is an explicit, iterative scheme which was developed 
by Matt Choptuik from the classic implicit Crank-Nicolson scheme jB], [7j. To solve equtation 
()2.2|) . we define the ICN-scheme as follows. First the iteration process is initiated: 

where S^{g^) = ^^'2 — — is the centered second order accurate finite difference approximation 
to the second order spatial derivative, g^ and are determined at mesh points Xi = iAx, 
t" = nAt and NLT {g^, Kf) is the non-linear term (i.e for eq. (Q NLT {g^, K^) - 
The scheme is then iterated: 



9n 



k\^^ = + At ( ^^^i 1 + NLT 




(/f^ =g^ + At 



n 



(2.15) 



j G [2,jmax]), and finally the dependent variables at the next time step are: 



jy-{n+l) jy-ijmax) 

i i 

(n+1) _ Umax) 

yi yi 



(2.16) 



As shown in and the optimal number of iterations, for the scalar wave equation is jmax = 
3^, which is the scheme that we will investigate in this paper. This scheme can be shown to be 
second order accurate in both time and space by a Taylor series expansion P,[21- 

^Note that different autliors count the number of iterations in different ways, some do not count the first step as an 
iteration and hence state that the optimal number of iterations is 2. 
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2.3.2 The 3'rd order Runge-Kutta scheme 

To solve eq. ()2.2|) . the 3'rd order Runge-Kutta scheme (RK3) is defined as follows [S]: 



K. =K,+ 



9^' = 97 + 



(1) , /I (2) , (3) 

9i + 4:9i ' +9i 



(2.17) 



where 



K, 



(1) 

i 



(2) 



9P 



At [5' {g^) + NLT{glKr 
At [Kn 



At 

At 

At 

At 



^2 ( „n I 9i 



(1) 



9- + 



+ NLT\g^ + ^,Kr + ^ 



K. 



(1) 



S'[g^- gf^ + 2^r^ )+NLT{g^- gl" + - + 2K 



.(2) 



(1) 



'(1) 



'(2) 



i^f - Kf'^ + 2K 



.(2) 



(2.18) 



This scheme can be shown to be second order accurate in space and third order accurate in time 
by a Taylor series expansion. 



2.3.3 The 4'th order Runge-Kutta scheme 

To solve eq. ()2.2|1 the 4'th order Runge-Kutta scheme (RK4) is defined as follows |H]: 



K. 



n+l 



6 



9i 9i ' g 



(2.19) 
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where 



K. 



(1) 



9i 



(1) 



K. 



(2) 



K. 



(3) 



,(3) 



K. 



(4) 



9i 



(4) 



At [6^g-) + NLTig:,Kr)] 
At [Kr] 



At 
At 
At 

At 

At 
At 



5' 



97 + 



9i 



(1)" 



K, 



(1)' 



5'\9- + 



9i 



(2)' 



+ NLT\gl 



g^^ 



(2) 



( ^^r + ^ ]+NLT{ g^ + i^r + ^ 



,(3) 



'(3) 



(3) 



(2.20) 



This scheme can be shown to be second order accurate in space and fourth order accurate in 
time by a Taylor series expansion. 



2.3.4 Crank-Friedrichs-Levy Nonlinear scheme (CFLN) 

Another approach to solving eq. (j2.2p numerically is to notice that without the non-linear term, 
equation ()2.2|) is a scalar wave equation which can be solved by the classic explicit scheme by 
Courant, Friedrichs and Levy l^(see also cp. 10 in Tj): 



-^^97 + 9. 



n-l 



At2 

We can cast this scheme into a first order form: 



^IVi - 2^r + gU 
Ax^ 



Kr- = Kr-^+At{5\g-)) 



9^' = 97 



97 + At< 



(2.21) 



(2.22) 



Now, in order to use this scheme as a basis for solving equation ()2.2|) . we must add a term to 
evolve the non-linear part, moreover, we wish to do this with second-order accuracy at the grid 
points (xj,t") to ensure overall second-order accuracy of the scheme. We do this by evaluating 
the non- linear term using the following predictor-corrector style approach 0: 

f>n+i 



Kr^+At{S\g^) + NLT{g^,K: 
At 

Y 



(2.23) 



9^ 



n+l 



97 + ^tKl 
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This scheme can be shown to be second order accurate in both space and time by a Taylor 
series expansion. An advantage of the scheme is that the second-order accuracy is achieved with 
a relatively small number of right hand side operations, to calculate one time step CFLN requires 
2 evaluations of the non-linear terms, compared to 3 evaluations required for the ICN and RK3 
schemes and 4 evaluations required for the RK4 scheme. It is also noted that the scheme is 

n-i 

staggered in time, i.e. initial values are required at points g'j_^ j ^. 



2.4 The Von Neumann analysis 

To investigate the properties of the numerical schemes we use a Von Neumann analysis following 
a standard approach [ZI,[01: AH the schemes can be represented by the following evolution 
operator: 

U^^' = ^Wli') (2-24) 

where Uj'i is a set of dynamical variables, n and j are temporal and spatial indices respectively 
and / enumerates the dynamic variables. In our case, the operator SJi may depend on all 
components of UJi^i, at grid points (j = j — 1, j, j + 1) corresponding to a time layer n (or a 
combination of n and n — ^ in the staggered case of CFLN scheme). 
Perturbing U^i as 

U^,i = U-,i + 8Ul,, (2.25) 

where f/j^^ is the background solution, substituting equation ()2.25j) into (j2.24j) and doing a Taylor- 
expansion, we obtain 



= E Wl^^^l^' + O {{^Ul,)') (2.26) 
Assuming a perturbation of the form 

= ^^e-^^^""^, (2.27) 

where < /c < is the perturbation wave number and I = V— 1, and substituting ()2.27p into 
flT^ we obtain 

^ E T^^I^e'^'"''^^'' = Gl^^i'^V (2.28) 

The amplification matrix G^i i, 

Qgn 

nn _ hi I(j-j')Axk /r, r,q\ 

holds information about properties of the numerical schemes. 

For a linear finite difference equation, the amplification matrix only depends on the dis- 
cretization parameters At, Ax and the wave number fc, and needs to be calculated only once in 
order to know properties of the scheme at all times and grid points. For non-linear schemes the 
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amplification matrix depends on f/"; and must in principle be evaluated at all spatial points and 
at all points in time. 

The generally complex eigenvalues of the amplification matrix, Aj, hold information about 
the amplification and speed of a given perturbation for a single time step which can be extracted 
by calculating the modulus and argument of the eigenvalues respectively. The classic condition 
for numerical stability is that the spectral radius (i.e. the largest modulus of the eigenvalues) 
of the amplification matrix is less than or equal to 1 for all wave numbers [7], however this 
excludes the possibility of growing exponential solutions. A less strict stability requirement (the 
Von Neumann stability criterion) which allows for exponentially growing solutions is that the 
eigenvalues must satisfy 

|A| =max|Ai| < 1 +0(Ai) (2.30) 

i 

for all /c |7j. In this paper, we will refer to the spectral radius |A| as the amplification factor. We 
presents amplification factors as a function of k ■ Ax which runs in the range k ■ Ax G [0, vr] with 
k ■ Ax = 71 corresponding to the Nyquist frequency. 

It should be noted that due to the local linearizion involved in calculating the eigenvalues for 
non-linear schemes, the Von Neumann analysis is only locally valid. This also means that the 
analysis cannot be trusted for non-local wave modes, i.e. small wave numbers. However, we are 
most interested in high wave numbers, as experience tells us that numerical instabilities usually 
arises first at the Nyquist frequencies. Hence for analyzing the stability of a scheme one searches 
for conditions under which equation ()2.30p is valid, focusing on the Nyquist frequency, which 
usually reduces to a restriction on the relationship between At and Ax, known as the Courant 
number a = 4f • 



3 The linear scalar wave equation 

Before studying the non-linear schemes, we briefly summarize stability, dissipation and dispersion 
properties of the schemes applied for a scalar linear wave equation 

d^g{x,t) ^ d^g{x,t) 

dt^ dx^ ^ ' ' 

Table El shows which Courant intervals satisfies the Von Neumann stability criterion ()2.30p 
for the Nyquist frequency as calculated by a Von Neumann analysis in the limit At — > 0. 



Scheme 


Stable Courant interval 


ICN 


< a < 1 


RK3 


< a < 


RK4 


0<a<V2 


CFLN 


< a <1 



Table 1: Courant intervals which satisfies the Von Neumann stability criterion (|2.3fl)) for the linear wave 
equation ()3.1() 



9 



1.01 



0.99 - 

0.98 - 

0.97 - 

0.96 - 

0.95 - 

0.94 - 

0.93 - 

0.92 - 

0.91 -IS" 
(3) 

„„ LSi- 

0.9 I 

0.1-31 0.2-ji 0.3'7i 0.4-7[ 0.5-71 0.6-31 0.7-31 0.8-31 0.9-3[ n 0.1-31 0.2-31 0.3-31 0.4-jt 0.5-3t 0.6-31 0.7-31 0.8-31 0.9-3[ n 

k-AX k-Ax 

(a) Courant number a — 0.50 (b) Courant number a — \/0.75 

Figure 1: Amplification factor as a function of k-Ax at Courant number a) a = 0.50 and b) a = VO.75, 
respectively, for the 4 schemes investigated. Legend is (1) = ICN, (2) = RK3, (3) = RK4 and (4) = 
CFLN. 

FiguresHshows the amplification factors as a function of k-Ax for Courant numbers a = 0.50 
and a = respectively, the latter corresponds to the maximal stable Courant number for 

the RK3 scheme (cf. table E)). These figures are representative of the dissipative behavior of the 
four schemes. The plot (and all plots in this section) is valid for all choices of Ax. 

For zero wave numbers all schemes are non-dissipative, but with increasing wave numbers 
and increasing Courant numbers the method of lines schemes shows monotonically increasing 
dissipation, with ICN being the most dissipative scheme followed by the RK3 and RK4 scheme 
respectively. For very high wave numbers the dissipation for the method of lines schemes shows 
a non-monotonic behavior as can be seen by comparing figures ^ We see that the ICN scheme 
is still, generally, the most dissipative scheme, but as Courant numbers are increased towards 
the stable limit, the dissipation at the Nyquist frequency vanishes and the maximal dissipation 
is seen for smaller wave numbers. Figure El shows the dissipative behavior of the schemes at 
the Nyquist frequency as a function of Courant number. This shows how the Nyquist frequency 
has maximal dissipation at around a ~ 0.8 ■ amax (where amax is the maximal stable Courant 
number for the corresponding scheme), after which the dissipation at the Nyquist frequency goes 
to zero just before numerical instability sets in. 

The CFLN scheme, in contrast, is completely non-dissipative for all wave numbers for all 
stable Courant numbers. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.f 
Courant number, a 



Figure 2: Largest eigenvalues at the Nyquist frequency (kAx = it) as function of Courant numbers for 
the 4 schemes investigated. Legend is (1) = ICN, (2) = RK3, (3) = RK4 and (4) = CFLN. 



Figures El compares the dispersion errors for the four schemes for various Courant numbers as 
calculated from an Von Neumann analysis for equation ()3.1|) . For small Courant numbers (fig. 
|3(a)| ), it is seen that the schemes behaves quite similarly (in fact, for a — > the dispersion errors 
for the four schemes converge to the same line), but a closer examination shows that the CFLN 
scheme has the smallest dispersion error, then follows the RK3 and RK4 schemes respectively 
and finally the ICN scheme which exhibits the largest dispersion errors. Also, dispersion errors 
are largest for large wave numbers and going to zero in the limit A; — > for all schemes. 
The dispersion errors for the ICN and RK4 schemes are increased monotonically for increasing 
Courant numbers at smaller wave numbers, while the dispersion errors for these schemes show 
some non-monotonic behavior at high wave numbers. As can be seen from figure |3(d)[ the ICN 
scheme in the limit of a = 1.00 exhibits positive dispersion errors for high wave numbers, i.e. 
the numerical solution is propagating faster then the analytic solution^. 

The dispersion errors for the RK3 and CFLN schemes conversely are minimized for their 
respective maximal stable Courant numbers. The CFLN scheme in this limit has zero dispersion 
errors, while the RK3 still has a non-vanishing (but minimized) dispersion error. 



^However, this is unlikely to be of numerical importance due the high damping for ICN in the high wave number / 
high Courant number limit. 
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k-Ax 
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k-Ax 



(a) Courant number a = 0.25 



(b) Courant number a = 0.50 
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0.1-Jr 0.2H 0.3!i 0.4jt 0.5-Jt 0.6-ji 0.7?i 0.8Tt 0.9-JE 7t 0.1-TZ 0.2-ji 0.3!i 0.4-jt 0.5-K O.Gtz 0.7-ji O.Sti 0.9jt jt 

k-AX k-Ax 



(c) Courant number a = 0.75 



(d) Courant number a — 1.00 



Figure 3: Wave speeds for the four schemes for various Courant numbers relative to the analytic wave 
speed. Legend is (1) = ICN, (2) = RK3, (3) = RK4 and (4) = CFLN. 



4 Von Neumann stability analysis of the non-linear wave equation 

We are interested in the behavior of the schemes in the non-hnear regime. In this section we 
present the results of a Von Neumann stabihty analysis of the four schemes applied to a local 
linearizion of the non-linear wave equation, equation ()2.2j) . presented in section EH) 
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k-Ax 



(a) Decaying solution 12.311 . 



(b) Growing solution 12.311 . 



Figure 4: Amplification factors for solution ()2.3|) as a function of k ■ Ax. (a) is decaying (C2 = —1) 
and (h) is growing (C2 = Ij- Courant number is a = 0.5 and Ax = j^. Legend is (1) = ICN, (2) = 
RK3, (3) = RK4 and (4) = CFLN. 

4.1 The spatial constant solutions 



Figures |3] shows a typical plot of amplification factors versus k ■ Ax for the spatially constant 
solution (j2.3p of equation (|2.2|) for decaying and growing solutions respectively, for the four 
schemes. Solution parameters are Ci = 2 for both plots and C2 = — 1 and C2 = 1 for figure |4(a)| 
and |4(b)| respectively. From the perturbation analysis in subsection 12.21 we expect to observe 
an amplification for all wave numbers when solution (j2.3p is decaying. For growing solutions 
we expect to see a more modest amplification at small wave numbers and damping at higher 
wave numbers. Looking at figures |3J we observe that at small wave numbers, all schemes agree 
with this prediction, i.e. we see a strong amplification at small wave numbers for the decaying 
solution and a smaller amplification for the growing solution. At high wave numbers the CFLN 
scheme agrees well with the perturbation analysis and we see an amplification at the Nyquist 
frequency, while the method of lines schemes all shows various degrees of damping, consistent 
with the results from section |21 The amount of damping is dependent on the Courant number 
as in the linear case, but also upon the spatial step size. Ax. By choosing a smaller Aa;, for 
a fixed Courant number. At decreases proportionally, by which any amplification also becomes 
smaller. However, with the amplification factor moving closer to |A| = 1, the damping for the 
method of lines schemes are shifted proportionally, i.e. the method of lines schemes may go from 
showing amplificative behavior to damping behavior at the Nyquist frequency, with the choice 
of a smaller Ax. We see this effect in figure |4(a)| for the ICN scheme. If we for the same solution 
had chosen a sufficiently small Ax, the other method of lines schemes would also have shown 
damping behavior at the Nyquist frequency. The qualitative behavior of the CFLN scheme, in 
contrast, is unaffected by the choice of Ax for reasonably small Ax. For very large Ax, all the 
schemes shows abnormal behavior which is caused by the fact that the base solution in this case 
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changes much more rapid then the temporal resolution allows for. 



4.2 The exponential solutions 

Figures El shows a typical plot of the amplification factors for the four schemes under investigation 
for the exponentially growing and decaying solution (j2.4j) to equation (j2.2j) as a function of k ■ Ax 
with solution parameter C = 1. We see that the plots are very similar to figures El For the 
decaying solution (figure 5(a) ) we see that CFLN is an amplificative scheme whereas the method 



of lines schemes displaying certain amount of dumping. We note that in order to reproduce a 
correct behavior of perturbed analytic solutions the scheme must have the amplification number 
greater than one. However, in numerical simulations it may be better to dump growing large 
wave number (short wavelength) perturbations instead of trying to faitfully reproduce them. 

For small wave numbers both the decaying and growing solutions are indicating amplificative 
behavior. We note that solution (j2.4j) is spatially dependent, hence the spatially constant per- 
turbation analysis from subsection 12.21 is not valid in this regime. Nevertheless, comparing with 
figures m where we see (and expect) the same behavior it seems plausible that the amplificative 
behavior seen in the figures is an expected behavior. 

As a closing remark to this section, we note that all schemes and solutions has been tested in 
the limit At ^ to investigate the stability properties of the schemes in the non-linear regime 
and we have found that all schemes agree with the stability intervals for the linear case (table 
121) for all solutions. 
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(a) Decaying solution 12.411 



(b) Growing solution 12.411 



Figure 5: Amplification factors for solution H2.4() as a function of k ■ Ax. (a) is decaying and (b) is 
growing. Courant number is a = 0.5, Ax = and C = 1. Legend is (1) = ICN, (2) = RK3, (3) = 
RK4 and (4) = CFLN. 
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5 Numerical Tests 



We have done numerical testings to verify the analytic results presented in the preceding sections. 
Figure ini shows the convergence of solution (j2.3|) to equation (|2.2|) . Plotted on the vertical axis 
is the absolute error of the central point in the domain between a simulation with the spatial 
step specified on the horizontal axis and a reference simulation. The reference simulation is a 
high resolution simulation with a resolution twice that of the leftmost point in the plot. The 
constants in solution ()2.3|) were set to Ci = 2 and C2 = —1, the simulations were run for a time 
t = 7r/4 in a domain of size 27r and the Courant number of the simulations was fixed at a = 0.50. 
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Figure 6: Absolute error between simulation at specified Ax vs. simulation at Ax = 1^ for Courant 
number a = 0.50. Legend is (1) = ICN, (2) = RK3, (3) = RK4 and (4) = CFLN. 



From the figure we see that the CFLN and ICN schemes are showing second order convergence, 
the RK3 scheme shows third order convergence. The RK4 scheme shows fourth order convergence 
at large Ax, but flattens out a around 10~^^ and for small step sizes. The flattening is caused 
by machine precision errors affecting the solutions. The convergence rates are what should be 
expected for spatially constant solutions. This solution does not introduce any truncation errors 
associated with spatial discretization. The figure thus shows convergence rates in agreement with 
the truncation error due to the temporal discretization of the schemes (second order CFLN and 
ICN, third order RK3, and fourth order RK4). Identical convergence tests for the exponential 
solutions ()2.4|) shows second order convergence for all schemes due to the second order spatial 
finite difference operator we have used in the schemes to calculate spatial derivatives. 

To test analytic predictions of the behavior of a perturbed solution we have made numerical 
tests of base solutions perturbed by small amplitude sinusoidal perturbations, g{x, t) = go{x, t) + 
g{x,t,k), with wave number k. Figure [7| is identical to figure El except that now base solution 
(j2.3|) is perturbed by a moving sinusoidal wave with wave number of A; = 16 (corresponding 
to the Nyquist frequency for the rightmost point) and initial amplitude Aq = 10~®, all other 
parameters are as for figure IHl 
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Figure 7: Absolute error between perturbed simulation at specified Ax vs. simulation at Ax = for 
Courant number a = 0.50. Perturbation is a moving sinusoidal wave with wave number k = 16 and 
initial amplitude = lO'^. Legend is (1) = ICN, (2) = RK3, (3) = RK4 and (4) = CFLN. 



We see from the figure that the perturbed solution is converging for all schemes. For higher 
resolutions the convergence is at least second order. Similar convergence rates are seen for 
analogous simulations with other base solutions. 

Figure IHl shows the ratio of the final amplitude (At) of a perturbation at time t = n/Ato its initial 
amplitude {Aq = 10"®) as a function of the resolution of the perturbed wave for simulations with 
base solution 12.31 with Ci = 10 and C2 = —1. Courant a = 0.50 and Ax = ^ was used in 
all simulations. The rightmost data point corresponds to a perturbation being at the Nyquist 
frequency. The analytical prediction is that 

A=./4M(0) = exp --^^1^— ^ 1.0435 

according to ()2.12j) . This prediction is also shown on the figure as line (0). 

From this figure we see that at high resolutions all schemes produce the same growth of amplitude 
of perturbations, which is in agreement with the analytical prediction. As the resolution decreases 
the results for all schemes for all schemes begin to deviate from the analytical value. Among 
the method of lines schemes, the RK4 is the most accurate and is able to reproduce the correct 
behavior of the perturbation up to the resolution of 16 grid points per wavelength. CFLN also 
requires only 16 grid points per wavelength to reproduce the correct behavior. The other two 
schemes requires from 32 to 64 grid points per wavelength. The major difference between the 
CFLN and the method of lines schemes is that at low resolutions, CFLN is amplifying the 
perturbations more than it should according to the analytical predictions, whereas the method 
of lines schemes are damping. The RK4 is the least damping and the most accurate of these 
schemes. The behavior of the schemes at low resolutions is fully consistent with the results of 
the Von Neumann stability analysis, presented in section 14.11 According to that analysis the 
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Figure 8: Ai/Aq of sinusoidal perturbation at t = 7r/4 for base solution (|2.|-i|) versus various wave 
numbers for a fixed Ax = ^ and a = 0.50. Legend is (0) = Analytic prediction, (1) = ICN, (2) = 
RK3, (3) = RK4 and (4) = CFLN. 



CFLN scheme must be amplifying at all wavelengths whereas the method of lines schemes must 
be damping. 

6 Discussion and conclusion 

In this paper we studied the properties of four numerical schemes, CFLN, ICN, RK3 and RK4, 
applied to a non-linear scalar wave equation ^2.1\ . This equation has a number of non-trivial 
analytic solutions whose properties, including stability, were studied and summarized in section 
El We carried out the Von Neumann stability analysis of the schemes and studied their phase 
and amplitude errors. Finally we carried out numerical experiments and compared the results of 
those experiments with the perturbation analysis of the equation and the Von Neumann stability 
analysis of the schemes. 

We find that all four schemes presented in this paper are stable and converge with the second 
order accuracy. The stability range for the schemes were determined and are presented in table 
121 Those ranges are valid for both linear and non-linear solutions. 

For non-linear schemes we checked that the amplification factor |A| is less than 1 + 0(At) for 
sufficiently small At. 

With respect to the dissipation errors, we find that the CFLN has the least amount of dissipation. 
The ICN scheme has the largest dissipation errors. The RK3 and RK4 are intermediate. For the 
method of lines schemes we find that damping is behaving non-monotonically with the increase 
of the Courant number. 

With respect to phase errors, the schemes can be arranged in the sequence CFLN RK3, RK4 
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and ICN, with the CFLN having the least amount of errors and ICN having the largest. 

We find that CFLN scheme requires the least amount of operations per time step, whereas the 
RK4 requires the largest amount. 

Ideally a numerical scheme is preferable which has the minimal phase and amplitude errors, 
and is computationally inexpensive. It is also of practical importance to have a scheme which 
will damp the high frequency perturbations at and close to the Nyquist frequency. Otherwise 
truncation errors at the highest frequency will remain within the computational domain, may 
be amplified and may eventually spoil the solution. None of the schemes discussed in this paper 
satisfy all those criteria. We think that the RK4 scheme should be preferred to other schemes 
because of its damping properties at the Nyquist frequency and the minimal amount of errors 
consistent with this property. On the other hand, the amount of dissipation in the RK4 scheme 
is dependent upon Courant number and is not controllable. If not for the damping properties, 
the CFLN scheme is the most accurate and cost effective. It would be of great interest to 
develop a version of this scheme with a controllable filter, for damping the minimum amount of 
perturbations at a narrow range of frequencies near the Nyquist frequency. 
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